home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / src / filter.cc < prev    next >
C/C++ Source or Header  |  1996-11-03  |  7KB  |  296 lines
  1. /*
  2.  
  3. Copyright (C) 1996 John W. Eaton
  4.  
  5. This file is part of Octave.
  6.  
  7. Octave is free software; you can redistribute it and/or modify it
  8. under the terms of the GNU General Public License as published by the
  9. Free Software Foundation; either version 2, or (at your option) any
  10. later version.
  11.  
  12. Octave is distributed in the hope that it will be useful, but WITHOUT
  13. ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  14. FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  15. for more details.
  16.  
  17. You should have received a copy of the GNU General Public License
  18. along with Octave; see the file COPYING.  If not, write to the Free
  19. Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
  20.  
  21. */
  22.  
  23. // Based on Tony Richardson's filter.m.
  24. //
  25. // Originally translated to C++ by KH (Kurt.Hornik@ci.tuwien.ac.at)
  26. // with help from Fritz Leisch and Andreas Weingessel on Oct 20, 1994.
  27. //
  28. // Rewritten to use templates to handle both real and complex cases by
  29. // jwe, Wed Nov  1 19:15:29 1995.
  30.  
  31. #ifdef HAVE_CONFIG_H
  32. #include <config.h>
  33. #endif
  34.  
  35. #include "defun-dld.h"
  36. #include "error.h"
  37. #include "oct-obj.h"
  38. #include "help.h"
  39.  
  40. extern MArray<double>
  41. filter (MArray<double>&, MArray<double>&, MArray<double>&);
  42.  
  43. extern MArray<Complex>
  44. filter (MArray<Complex>&, MArray<Complex>&, MArray<Complex>&);
  45.  
  46. template <class T>
  47. MArray<T>
  48. filter (MArray<T>& b, MArray<T>& a, MArray<T>& x, MArray<T>& si)
  49. {
  50.   MArray<T> y;
  51.  
  52.   int a_len  = a.length ();
  53.   int b_len  = b.length ();
  54.   int x_len  = x.length ();
  55.  
  56.   int si_len = si.length ();
  57.  
  58.   int ab_len = a_len > b_len ? a_len : b_len;
  59.  
  60.   b.resize (ab_len, 0.0);
  61.  
  62.   if (si.length () != ab_len - 1)
  63.     {
  64.       error ("filter: si must be a vector of length max (length (a), length (b)) - 1");
  65.       return y;
  66.     }
  67.  
  68.   T norm = a (0);
  69.  
  70.   if (norm == 0.0)
  71.     {
  72.       error ("filter: the first element of a must be non-zero");
  73.       return y;
  74.     }
  75.  
  76.   y.resize (x_len, 0.0);
  77.  
  78.   if (norm != 1.0)
  79.     b = b / norm;
  80.  
  81.   if (a_len > 1)
  82.     {
  83.       a.resize (ab_len, 0.0);
  84.  
  85.       if (norm != 1.0)
  86.     a = a / norm;
  87.  
  88.       for (int i = 0; i < x_len; i++)
  89.     {
  90.       y (i) = si (0) + b (0) * x (i);
  91.  
  92.       if (si_len > 1)
  93.         {
  94.           for (int j = 0; j < si_len - 1; j++)
  95.         si (j) = si (j+1) - a (j+1) * y (i)
  96.           + b (j+1) * x (i);
  97.  
  98.           si (si_len-1) = b (si_len) * x (i)
  99.         - a (si_len) * y (i);
  100.         }
  101.       else
  102.         si (0) = b (si_len) * x (i)
  103.           - a (si_len) * y (i);
  104.     }
  105.     }
  106.   else if (si_len > 0)
  107.     {
  108.       for (int i = 0; i < x_len; i++)
  109.     {
  110.       y (i) = si (0) + b (0) * x (i);
  111.  
  112.       if (si_len > 1)
  113.         {
  114.           for (int j = 0; j < si_len - 1; j++)
  115.         si (j) = si (j+1) + b (j+1) * x (i);
  116.  
  117.           si (si_len-1) = b (si_len) * x (i);
  118.         }
  119.       else
  120.         si (0) = b (1) * x (i);
  121.     }
  122.     }
  123.   else
  124.     y = b (0) * x;
  125.  
  126.   return y;
  127. }
  128.  
  129. extern MArray<double>
  130. filter (MArray<double>&, MArray<double>&, MArray<double>&,
  131.     MArray<double>&);
  132.  
  133. extern MArray<Complex>
  134. filter (MArray<Complex>&, MArray<Complex>&, MArray<Complex>&,
  135.     MArray<Complex>&);
  136.  
  137. template <class T>
  138. MArray<T>
  139. filter (MArray<T>& b, MArray<T>& a, MArray<T>& x)
  140. {
  141.   int a_len = a.length ();
  142.   int b_len = b.length ();
  143.  
  144.   int si_len = (a_len > b_len ? a_len : b_len) - 1;
  145.  
  146.   MArray<T> si (si_len, T (0.0));
  147.  
  148.   return filter (b, a, x, si);
  149. }
  150.  
  151. DEFUN_DLD (filter, args, ,
  152.   "usage: [y [, sf]] = filter (b, a, x [, si])\n\
  153. \n\
  154. y = filter (b, a, x) returns the solution to the following linear,\n\
  155. time-invariant difference equation:\n\
  156. \n\
  157.   a[1] y[n] + ... + a[la] y[n-la+1] = b[1] x[n] + ... + b[lb] x[n-lb+1],\n\
  158. where la = length (a) and lb = length (b).\n\
  159. \n\
  160. [y, sf] = filter (b, a, x, si) sets the initial state of the system, si,\n\
  161. and returns the final state, sf.  The state vector is a column vector\n\
  162. whose length is equal to the length of the longest coefficient vector\n\
  163. minus one.  If si is not set, the initial state vector is set to all\n\
  164. zeros.\n\
  165. \n\
  166. The particular algorithm employed is known as a transposed Direct Form II\n\
  167. implementation.")
  168. {
  169.   octave_value_list retval;
  170.  
  171.   int nargin  = args.length ();
  172.  
  173.   if (nargin < 3 || nargin > 4)
  174.     {
  175.       print_usage ("filter");
  176.       return retval;
  177.     }
  178.  
  179.   const char *errmsg = "filter: arguments must be vectors";
  180.  
  181.   int x_is_vector = (args(2).rows () == 1 || args(2).columns () == 1);
  182.  
  183.   int si_is_vector = (nargin == 4
  184.               && (args(3).rows () == 1 || args(3).columns () == 1));
  185.  
  186.   if (args(0).is_complex_type ()
  187.       || args(1).is_complex_type ()
  188.       || args(2).is_complex_type ()
  189.       || (nargin == 4 && args(3).is_complex_type ()))
  190.     {
  191.       ComplexColumnVector b = args(0).complex_vector_value ();
  192.       ComplexColumnVector a = args(1).complex_vector_value ();
  193.       ComplexColumnVector x = args(2).complex_vector_value ();
  194.  
  195.       if (! error_state)
  196.     {
  197.       if (nargin == 3)
  198.         {
  199.           ComplexColumnVector y (filter (b, a, x));
  200.  
  201.           if (x_is_vector)
  202.         retval (0) = octave_value (y, (args(2).columns () == 1));
  203.           else
  204.         retval (0) = y;
  205.         }
  206.       else
  207.         {
  208.           ComplexColumnVector si = args(3).complex_vector_value ();
  209.  
  210.           if (! error_state)
  211.         {
  212.           ComplexColumnVector y (filter (b, a, x, si));
  213.  
  214.           if (si_is_vector)
  215.             retval (1) = octave_value (si, (args(3).columns () == 1));
  216.           else
  217.             retval (1) = si;
  218.  
  219.           if (x_is_vector)
  220.             retval (0) = octave_value (y, (args(2).columns () == 1));
  221.           else
  222.             retval (0) = y;
  223.         }
  224.           else
  225.         error (errmsg);
  226.         }
  227.     }
  228.       else
  229.     error (errmsg);
  230.     }
  231.   else
  232.     {
  233.       ColumnVector b = args(0).vector_value ();
  234.       ColumnVector a = args(1).vector_value ();
  235.       ColumnVector x = args(2).vector_value ();
  236.  
  237.       if (! error_state)
  238.     {
  239.       if (nargin == 3)
  240.         {
  241.           ColumnVector y (filter (b, a, x));
  242.  
  243.           if (x_is_vector)
  244.         retval (0) = octave_value (y, (args(2).columns () == 1));
  245.           else
  246.         retval (0) = y;
  247.         }
  248.       else
  249.         {
  250.           ColumnVector si = args(3).vector_value ();
  251.  
  252.           if (! error_state)
  253.         {
  254.           ColumnVector y (filter (b, a, x, si));
  255.  
  256.           if (si_is_vector)
  257.             retval (1) = octave_value (si, (args(3).columns () == 1));
  258.           else
  259.             retval (1) = si;
  260.  
  261.           if (x_is_vector)
  262.             retval (0) = octave_value (y, (args(2).columns () == 1));
  263.           else
  264.             retval (0) = y;
  265.         }
  266.           else
  267.         error (errmsg);
  268.         }
  269.     }
  270.       else
  271.     error (errmsg);
  272.     }
  273.  
  274.   return retval;
  275. }
  276.  
  277. template MArray<double>
  278. filter (MArray<double>&, MArray<double>&, MArray<double>&,
  279.     MArray<double>&);
  280.  
  281. template MArray<double>
  282. filter (MArray<double>&, MArray<double>&, MArray<double>&);
  283.  
  284. template MArray<Complex>
  285. filter (MArray<Complex>&, MArray<Complex>&, MArray<Complex>&,
  286.     MArray<Complex>&);
  287.  
  288. template MArray<Complex>
  289. filter (MArray<Complex>&, MArray<Complex>&, MArray <Complex>&);
  290.  
  291. /*
  292. ;;; Local Variables: ***
  293. ;;; mode: C++ ***
  294. ;;; End: ***
  295. */  
  296.